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The high temperature many-body density matrix is fundamental to path integral computation. 
The pair approximation, where the interaction part is written as a product of pair density matrices, 
is commonly used and is accurate to order r^, where r is the step size in the imaginary time. Here 
we present a method for systems with Coulomb interactions in periodic boundary conditions that 
consistently treats the all interactions with the same level of accuracy. It shown that this leads 
^ ' to a more accurate high temperature solution of the Bloch equation. The method is applied to 

many-body simulation and tests for the isolated hydrogen atom and molecule are presented. 
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I. INTRODUCTION 



Quantum Monte Carlo (QMC) methods are frequently used to study interacting many-body systems in different 
^ fields of physics and chemistry when a high degree of accuracy is needed [l|. The description of correlation effects 
combined with favorable scaling properties of or better (N is the number of particles) make these techniques 
effective in many applications. Path integral Monte Carlo (PIMC) is unique among other QMC methods because it 
can describe quantum systems at finite temperature 0, 0, 13; S • The method is based on the thermal density matrix 
that characterizes the properties of a system in thermal equilibrium. 

Many applications of the PIMC method require an accurate treatment of Coulomb interactions in periodic boundary 
conditions (PBC). This includes all electronic structure simulations that describe electrons as individual particles. 
However, there is still no universally accepted method to compute the Coulomb propagator in PBC, which has led to 
unnecessary approximations resulting in less efficient and less accurate many-body simulations. 
fH . In this article, we describe an accurate approach to compute the Coulomb pair density matrix in a periodic system. 
Q ■ It can easily be generalized to other long range interactions [7|. The density matrix of a bosonic (B) or fermionic (F) 
O . system at temperature T can expressed in terms of an imaginary-time path integral. 



PB/F(R,R';/?) = ^^(±l)''y" rfR» p(R,Ri;r)p(Ri,R2;r) . . .p(RM-i,PR';r) , (1) 

00 

' where r is the time step r = (3/M with (3 = and fcs is Boltzmann's constant. particles in real-space 

1 representation, R — {ri, . . . , r^v}, are described in the canonical ensemble. Atomic units of Bohr radii and Hartree 
will be used throughout this work. 

Instead of requiring the propagator at temperature, T, path integrals rely a density matrix at much higher tem- 
perature, M X T. At high temperature, the many-body density matrix can be computed with good accuracy because 
I exchange effects as well as three-body correlations are negligible in the limit of high temperature. A novel method 
of constructing the high temperature density matrix (HTDM) for systems with long-range interactions with periodic 
boundary conditions is the focus of this study. 
Different approximations for the Coulomb HTDM have been advanced that all become exact in the limit of r 0. 
^ ' However, all PIMC simulations are performed at finite t. Therefore an accurate representation of the HTDM is 
very important. Its accuracy determines the maximum time step t one can use. A larger time step allows one to 
significantly cut down on the number of time slices in the path integral, M, and therefore improve the efficiency of 
many-body simulations. This gain in efficiency may be crucial in practical applications to perform accurate simulations 
at low temperature or for large systems. 

The HTDM can be used to define the potential action, U (R, R'; t). For a system of N distinguishable particles in 
real space representation this reads, 

JV 

p(R,R';T) = exp{-C/(R,R';T)} []po(r„r^;r) . (2) 

i=l 

Pq is the free particle density matrix in D dimensions, 

Po(r, r'; r) = (47rAr)-^/2 exp{-(r - r') V4At} . (3) 
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A is the a mass dependent parameter, fi^ /2m. U is commonly approximated as the product over the nonideal parts 
of all pair density matrices p(r,y, rj^ ; r), 

exp{-[/(R,R^r)}«exp|-^u(r,,,4,r)l =nf^^4^ ■ (4) 

This is called the pair approximation, and is accurate to order 0j S 0| ■ One essentially starts with a exact solution 
of the two-body problem when performing many-body simulations. This mean only 1 time slice is needed to study the 
hydrogen atom at any temperature and equilibrium ionization state but more are needed for the hydrogen molecule 
when three-body correlations are important. Exchange effects have been neglected in Eq.|4l which justified for a small 
time step r. In the full path integral, they enter as a sum of over permutations in Eq. [T] 

For a Coulombic system with periodic boundary conditions, these pair density matrices, p(r,r';T), are solutions of 
the two particle Bloch equation with the Ewald potential Q, 

^ = -Hp = p - VEw{r)p , (5) 

with the initial condition p(r,r';T = 0) = 5{r — r'). The reduced mass fiij = niimj / {nii + nij) enters through 
Ay = h^/2fj,ij. The pair density matrix can also be derived from the Feynman-Kac relation, 

p(r,r';r) =po(r,r';r) (e" ''(■"^*»)_^, , (6) 

where the average is to be taken over all free-particle (Brownian) paths from r to r'. For potential without negative 
singularities, this expression can easily be evaluated numerically for specific pairs of r and r' in order to verify the 
approximations to be discussed below. However, the results will always have a statistical uncertainty due to the finite 
sample of Brownian paths. 

The resulting Ewald pair action, MEw(r, r'; t) = ^ lii[jOEw/po)]7 determines the weight of the paths in Eq. [1] where 
r and r' represent the separations of pairs of particles i and j at two adjacent time slices, 

r^r,{t)-rj{t) , r' = r,(i + r) - r,(i + r) . (7) 

For a spherically symmetric potential, the pair density matrix depends on r and three spatial variables: the initial 
and final pair separations |r| and |r'| as well as the angle between them 9. Alternatively, it can be expressed in terms 
of the variables q, s and z, 

g^i(lrKlr'l), .^|r-r'|, z^|r|-|r'|. (8) 

For Coulomb potential, the dependence on z drops out The Ewald potential, however, has the symmetry of the 
periodic simulation cell which requires both the initial and final pair separation to be specified with respect to the 
cell. This implies that the pair density matrix for the Ewald potential depends on six spatial variables. This makes 
computation and storage of the corresponding action extremely awkward. 

Previous methods [13, |Tl| to deal with this difficulty have involved a break-up of the Ewald potential into a 
spherically symmetric short-range piece and a long-range remainder, 

VEwir)^V,,,X\v\) + V,,,Xr) ■ (9) 

The short-range piece has been treated numerically using the matrix squaring technique developed by Storer . In 
principle, it allows one to derive the exact action for a spherically symmetric potential but in practice the accuracy 
is controlled by numerical accuracy of the integration. Matrix squaring is performed on a grid and controlling the 
associated grid errors requires significant care [l3|. To treat the cusp condition at the origin accurately, the short- 
range part must include the singular part of the potential. Most simply, one can use the direct 1/r interaction term 
as short-range part. Alternatively, one can employ the optimized Ewald break-up method described in [3], which 
allows one to construct a short-range piece that always decays within the boundaries of the simulation cell. A detailed 
review of the break-up method and its accuracy is given in ^13.] . The primitive approximation, 



up.A.{r,r';r)^^[V{r) + V{r')] , 



(10) 
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provides a simple straightforward way to add the fong-range remainder, 

u{r,r';T) ^ u,.,.{r,r';T) + + Vi.r.{r')] ■ (11) 

but the random phase approximation has also been applied to the long-range action When we later refer to 

Eqs. [9landfTT| we assume that the l/r pot ential has been used as short-range piece. 

The break-up of the Ewald potential [id . [llj introduces an additional approximation to path integral computation. 
While it has been successfully in many applications [iM [3 [3 there is a need for improvement. The break-up is 
ad hoc and introduces some arbitrariness to the construction of the Ewald action. The accuracy remains high if cither 
potential piece is sufhciently smooth on the scale of the thermal de Broglie wave length, Xd = \/4tt\t. However, the 
convergence to the correct answer as function of cell size and temperature is difficult to assess. For these reasons, 
we present a method here that consistently treats the direct Coulomb interaction within the simulation cell and that 
with periodic images with the same level of accuracy. Our method thereby avoids introducing a short and long-range 
potential. The necessary pair density matrices are derived for the Coulomb potential, which makes it possible to 
utilize the large amount of analytic and numerical work available for this potential. The development of the method 
proceeds by three steps which are now described. 



II. METHOD 



A. Computation of the pair density matrix for the Coulomb potential 

It is first necessary to compute the pair density matrix for an isolated pair of particles. This can either be done 
with matrix squaring or from the sum over eigenstates. In the matrix squaring technique, one can take advantage of 
the fact that only s-wave contributions enter in case of the Coulomb potential as done by Storer [l^ . In the squaring 
technique, one starts from a high temperature expansion and then numerically squares the density matrix using. 



p{r, r'; 2r) = J dr" p(r, r"; r) p(r", r'; r) , (12) 



until the lowest temperature has been reached. This typically requires 10-30 iterations. Since matrix squaring is 
done on a grid in r and r' a grid error is introduced. For the Coulomb potential this usually requires a high number 
grid points near the singularity at the origin. This singularity also requires care when the density matrix is initialized 
using a high temperature expansion at the beginning of the squaring procedure. Despite these shortcomings, the 
squaring technique has the advantage that it can easily be applied to arbitrary spherical potentials [lOj . Alternatively, 
Schmidt and Lee developed an approach where kinetic and potential operators are applied repeatedly using Fourier 
transforms [l^]. This method has been adopted by J. Shumway to study a variety of systems with PIMC [2]| . 
Vieillefosse used an analytic power-series expansion method to construct the action [22l | . 

For potentials that do not exhibit negative singularities, the pair density matrix can also be derived from the 
Feynman-Kac formula, Eq. [6l This approach was applied to the Yukawa interaction in Ref . Q . It avoids introducing 
grid errors and readily yields diagonal matrix elements as well as the first term in an expansion for off-diagonal 
elements but it impractical for constructing a table the includes arbitrary off-diagonal elements. Matrix squaring 
would then be more appropriate [isj . 

In this paper, we rely on another approach and compute the pair density matrix by summation over eigenstates [2^ . 

p(r,r';r)=^e-"=Vs(r)^:(r'). (13) 

s 

It also avoids grid errors and can provide off-diagonal elements efficiently but requires that the eigenstates are known 
with high precision, which is the case for the Coulomb potential. Although the number of states to be considered 
increases with temperature but the summation can be performed accurately for the temperatures of interest as was 
shown in Ref. (23|] . 

For illustration. Fig. [T] shows examples of the electron pair potential action for the on diagonal case (r' = r) and 
the exchange case (r' = — r) over two decades of r. The cusp condition leads to a linear behavior near the origin. For 
large r, the diagonal action converges to the primitive action. The exchange term always has a higher action, which 
can be understood from the Feynman-Kac formula in Eq. [S] where one averages over all Brownian paths. While in 
the exchange case paths beginning at r must diffuse around the repulsive l/r potential to reach — r, in the diagonal 
case can avoid region of high potential more easily. 
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FIG. 1: Potential action divided by the inverse temperature t at t — 0.1, 1.0, 10.0 for an isolated pair of electrons as a 
function of x' with r — (|x'|,0, 0) and r' = (a;', 0,0). With increasing separation, the diagonal action (r = r') converges 
rapidly to the primitive action, r/r. The action in exchange case (r — — r') is consistently higher than the diagonal action, 
ii(r, — r, r) > u(r, r, r). 
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FIG. 2: The potential pair action u{q, s), shown for a proton-electron pair at t = 1 shown as function of q — i[|r| + |r'|]. The 
left graphs shows that, for diagonal configurations r = r', the error in the primitive approximation, —r/q, decays like q~*'^ ■ 
The right graphs shows that the off-diagonal terms decays like q '''^ for fixed s = |r — r'| = 2-\/Ar — \/2 (thick dashed line). 
The dot-dashed line shows the error in the primitive approximation for parallel r and r' (^ = 0). The solid lines shows this 
error in the case of |r| = |r'| and angle 9 > 0. 



To increase the efficiency of PIMC simulations, the HTDM is derived beforehand and tabulated. For a given time 
step T, the Coulomb pair density matrix is a function of the two variable q and s, u{q, s) = u(r,r';T). In order to 
reduce the storage, we fit the off-diagonal terms using an expansion of powers of (compare with ^]), 

n 

u(g,s)«w(<z,0)+^A,(Q)s2^, (14) 
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For each q, a least squares fit can represent the off-diagonal term with sufficient accuracy. The results for n — 1 are 
shown in Fig. [3l The expansion works well in PIMC because the kinetic energy prevents adjacent point on the path 
to exceed separation much larger than the thermal de Broglie wave length. The storage of complete set of off-diagonal 
terms is therefore not needed. For the Coulomb potential, the off-diagonal terms decay approximately like q^'^'^, 
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FIG. 3: Off-diagonal expansion coefficients for proton-electron pair action Ai and its r derivative A'l at r = 1. 

which is demonstrated by three curves in Fig. [2l The first compares the diagonal (s = 0) and off-diagonal action 
(|r| = \r'\ ^ q, s ^ 2\/At = \/2) at constant q. The magnitude of s represents a typical separation of adjacent points 
along the path which is mostly determined by the free particle action. The second curve in Fig. [2] shows the error 
in the primitive approximation for the same q and s parameters. The third curve shows the error in the primitive 
approximation for parallel vectors {9 = 0) with the exact action. The graph underlines that the sum over off-diagonal 
contributions from the image charge converges. Also for large separations, the primitive approximation works well, 
and off-diagonal contributions may be neglected. 

These properties is used in PIMC simulations, where one typically computes the off-diagonal contributions only 
for particle pairs within the simulation cell using the minimum imagine convention. To date, no efficient way to sum 
up all off-diagonal contributions has been proposed. However, if necessary, a summation over the required number of 
images can easily be included. 



B. Pair approximation for the image charges 

The pair approximation Eq. 3] is commonly used to approximate the action of a path. This expression follows from 
the Feynman-Kac formula Eq. [6l here written for a finite system of N particles. 



exp{-;7(R,R';T)} = /expi- /"dtV 
\ [ ^0 r<J 



QiQj 



> (15) 

R^R' 

{n-{-i^'Kl}>_, 
n gi})_„. 

J|exp{-M(ry,r^;T)} (18) 



exp 



Y,u{v,j,v[^-t)) (19) 
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The pair approximation enters in Eq. [T7] where it is assumed that interactions between pairs of particles may be 
averaged independently from the location of paths of the remaining particles. This approximation becomes increasingly 
accurate for large separations. It is exact at high temperature and, as stated before, the error is of 0{t^). 

Now we will describe how the pair approximation can be applied to the systems with periodic boundary conditions. 
In principle, one only has to add sum over the interaction with of all periodic images using the pair approximation, 
which assumes that the diffusion of paths of the image particles may also be averaged independently. Since the error 
in the pair approximation also decreases with separation, applying it to the periodic images is less of an approximation 
than one is using already by applying the pair approximation to the direct interaction within the simulation cell. 

However, Eq. 1191 converges conditionally since the Coulomb potential is long-ranged. This problem has been solved 
by introducing the Ewald potential which assumes charge neutrality guaranteed by a neutralizing background. Starting 
from the Ewald potential, we will now give an approximate expression for the corresponding Ewald action. Since 
the exact action converges to the primitive approximation for large separations (Fig. [2]), the quantum correction is a 
short-ranged function, 

Au{r, r'; r) = u(r, r'; r) - ^Q,Q, + ^) • (20) 

Since it decays faster than r^^ (Fig. [2]), the sum over all images converges. It should be noted that u(r,r';r) does 
not need to go to zero with the simulation cell. This allows us to construct an Ewald action using the primitive 
approximation and a sum of quantum corrections Au, 



L V* ^^^^^ ^^^^^^^^^^^ ^^^^^^^^^^^^ . .v^, 

UEw(r, r'; r) w ^QiQj[VEw(r) + VEw(r')] + ^"("^ + L, r' + L; • 

L 



-)+UBG , (21) 



which is the central result of this article. We will discuss diagonal matrix elements first, off-diagonal elements are 
derived in the next section. Note that it is periodic as well as symmetric in the arguments r and r', and satisfies the 
cusp condition, i.e. the 1/r singularity in the Coulomb potential as r ^ is canceled in the Bloch equation by a term 
arising from the kinetic energy operator. 

The existence of the background term ubg can be understood in two ways. First, by introducing the residual of 
Bloch equation. 



i?(r,r';r)^ ^ 



p(r,r';r) 



^ - Ay Wl + V 

OT 



p(r,r';T). (22) 



The residual derived from the primitive action leads to the following simple expression on the diagonal, 

i?(r,r;r) = ^V?y-^(Vr^)^ . (23) 

For the Ewald potential, the presence of a neutralizing charge background (V^Vew = 47r/r2) gives rise a constant 
that is independent of r. An additional correction for it can be derived assuming a locally harmonic potential 3] at 
the point r* = (L/2, L/2, L/2), 

^^BG = ^Q.Q,^-5]Aii(r*+L,r*-|-L;r) . (24) 

The first term is of order O(t^) and is therefore not part of the primitive approximation. A constant shift in the 
action does not affect PIMC sampling, however, the r derivative of ubg enters in the estimator for kinetic energy. The 
second way to derive the background term is by numerically calculations using the Feynman-Kac formula[6l In Fig. HI 
the difference of the various approximations for the Ewald action are compared with the computed exact values. The 
graphs shows that the error in the primitive approximation approaches a constant at large r while it decays to zero for 
a pair of isolated particles. The difference primarily arises from the presence of the neutralizing background, which 
was corrected by introducing ubg- However, it should be stressed that, t = 1 in Fig. [4] is much larger than typical 
time steps used in simulations with electrons and nuclei. In such simulations, the background contributions are less 
importance because they scale like r^. 

In analogy to Ewald approach (Eg. IA2p . the action of a many-body system can be written as, 

f/(R, R'; ^) = E "EW., {r^J , ; t) -I- ^ t Q^^m ■ (25) 

i>j i 

where the charge factors have been included in the action terms. The efficient evaluation of uew is discussed in 
Appendix [X] 
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FIG. 4: Various approximations for the Ewald action for pair of electrons at r = 1 are compared on the diagonal as a function 
of |r| with the exact action u, which was either computed analytically or using the Feynman-Kac formula. The left graph 
shows action difference along the (1,0,0) direction in a cubic cell of 20 a.u. size. The right graph repeats all those curves with 
open symbols for a smaller cell of 5 a.u. but results for the (1,1,1) direction have been added with filled symbols. The error 
in the primitive approximation up. a. does not decay to zero, underlining contributions of the neutralizing charge background 
(dot-dashed line). 



III. RESULTS 



Tables U and |TT] give the density matrices and their r derivative for isolated proton-electron and electron-electron 
pairs at a time step of 0.125 Hartrees"^, which is typical for a PIMC simulation. In the notation of Eq. [TH the 
diagonal and first order off-diagonal expansion coefficients are tabulated. 

As a first test of the isolated pair action, a single hydrogen atom at a temperature T = 0.025. This temperature 
was chosen low enough for the atom to be in the ground state since the relative occupation probability of the first 
excited state is 5 x 10~^. Path integral Monte Carlo estimates of the potential energy agrees with the exact value of 
— 1 within the statistical error bar of the Monte Carlo integration, here 10"''. While the evaluation of the potential 
energy involves only the action, the total energy estimator requires also the t derivative. The computed total energy 
agreed with the exact value of —1/2 to a relative accuracy of 7 x lO"**. Both potential and kinetic energy estimators 
are necessary to calculate the pressure, P, from the virial theorem, 3 (P) fl = 2 (K) + (V) . 

As a second test for the accuracy of the computed pair density matrix, the kinetic and potential energy for an 
isolated hydrogen molecule was computed as a function of proton separation R. Fig. [5] shows very good agreement 
with the exact groundstate calculation [23|. The molecular binding energy of 1.174448 as well as the location of the 
minimum at i? = 1.4008 are well reproduced (see inset of figure[5]). The virial theorem for a diatomic molecule [24 12^ 
reads, 

2K + V = ~'R^. (26) 
dR 

It goes to zero at the equilibrium bond length as well as in the limit of i? ^ oo as shown in the figure. 

A more detailed analysis of the comparison with the groundstate calculation for R fixed at the equilibrium bond 
length shown in Fig. [5] reveals that the finite time step in the path integral determines the accuracy of the calculation. 
The temperature T — 0.025 was chosen low enough to compare with the groundstate result as the check with 
T = 0.0125 shows. The internal energy converges faster to the exact groundstate result than the virial term 2K + V. 
For T < 0.25, we find that the energy deviates less than 10~^ from the exact groundstate energy. For r = 0.0625, the 
energy deviated by (2.4 ± 1.7) x lO^"'. For this time step, we obtained 2K + V = (9± 7) x lO"**, which corresponds 
to an residual inaccuracy in the pressure equivalent to the pressure of an ideal molecular gas at T = 90 ± 70 K, which 
is a significant improvement in accuracy over [IZ] and more than sufficient for the majority of PIMC simulations of 
hot, dense hydrogen [16, 26]. This concludes the accuracy analysis for isolated systems of particles. 

The accuracy of the constructed pair density matrix in a periodic system is analyzed by calculating the error in 
Bloch equation [5] for the Ewald potential. The optimized Ewald break-up of Eqn. IA7I for a cubic simulation cell of 
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FIG. 5: Virial term 2K + V (kinetic energy K and potential energy V) and total internal energy E are shown for an isolated 
hydrogen molecule, for which the nuclei have been kept fixed at separation R. The PIMC results have been obtained from 
simulations at T = 0.0125 with 320 time slices. The solid lines are the exact groundstate results of [23 ]. 
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FIG. 6: Total internal energy E and virial term 2K -\-V for an isolated hydrogen molecule are shown for different temperatures 
as a function of the time step in the path integral. The nuclei have been kept fixed at the equilibrium separation of 7? = 1.4008. 



size L = 5 using 20 /c-shells is given in Tables IIIII and IIVI 

Based on the residual in Eq. 1221 we define the accuracy parameter, 

I{L-t)^ jj dvdv'p{v,v'-T)\R{v,v'-T)\l jjdvdv'p{v,v'-T), (27) 

that characterizes the accuracy of the particular approximation to the density matrix. Fig. [7] shows this accuracy 
parameter for the periodic proton-electron pair density matrix as a function of cell size L. This present method of 
treating the image charges in the pair approximation is contrasted with treating them in the primitive approximation, 
Eqn. [Til with the long range term Vi,r.{v) = Vew{'^) — ^/f- In the limit of large L, both methods are in agreement 
and show a rapid decay of the residual indicating that the Bloch equation is satisfied with increasing accuracy. For 
small L, the pair approximation yields a smaller residual demonstrating that the method proposed here leads to a 
more accurate high temperature density matrix for periodic systems. 
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We have applied Eq. [51] to PIMC simulations of hydrogen, helium, hydrogen-helium mixtures [13, [H, [H, [13] , and 
to study quantum effects in the one-component plasma [tI, 'sT', 's^l . Fig. [8] demonstrates the quality of the construct 
density matrices in PIMC simulation of a strong coupled Coulomb system. Despite significant quantum properties 
the particles arrange in a Wigner crystal. The presented time step analysis shows that a time step of 400 or less is 
needed to obtain converged structural properties such as pair correlation functions. This serves as a benchmark for 
future methods to construct density matrices for PIMC that would then need to demonstrate that the same level of 
accuracy can be achieved with larger time steps. 




L (a.u.) 

FIG. 7: Residual of the Bloch equation I{L;t) in periodic boundary conditions is shown as function of cell size L for the 
proton-electron pair density matrix at two different temperatures (open symbols — 8 and full symbols = 0.5). The 
pair approximation for periodic images (circle) is compared with treating the images in the primitive approximation (squares). 




FIG. 8; The pair correlation function g{r) for the homogeneous electron gas in state of a Wigner crystal [3, lisl . Isil . [3^ for 
coupling parameter F = 200 and quantumness parameter ?7 = 1 (/3 = 40000, = 200) are shown from simulation with different 
time steps t. The right graph shows how the maximum value of the g{r) converges as a function of r. Good convergence is 
reached for r = 400 requiring simulations with 100 time slices. 
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IV. CONCLUSIONS 



The method proposed here to compute the Coulomb pair density matrix in periodic boundary conditions has two 
advantages over previous approaches. It starts from the more accurate Coulomb pair density matrix, which reduces 
the numerical errors in the short-range behavior of the computed pair density matrix and circumvents a break-up 
of the potential. Although it does not derive the pair density matrix for the Ewald potential exactly, it consistently 
uses the pair approximation for the action which is an intrinsic assumption in most many-body path integration 
simulations. 

In addition to describing this method, quantitative test results for the isolated hydrogen atom and molecule as 
well as a test for the periodic system based on the residual error in the Bloch equation were presented to assess the 
accuracy of the high temperature density matrix. It is suggested that these tests should form a basis for comparing 
any proposed alternative methods. 



APPENDIX A: EWALD POTENTIAL 



The Ewald method to describe the potential of charged particles in periodic boundary conditions can be adapted 
in the following way to calculate the action. The potential energy of a system of N charges Qi interacting via the 
Coulomb potential K(r) = l/|r| is given by [s^ . 

^ = E E Q'Q^ ^-(^^^- + L) + ^ E E ^-(i^) ■ (Ai) 

i>j L Lt^O i 

where = — Vj and L is a lattice vector. This expression converges conditionally for charge neutral systems 
Qi = 0)- Under the assumption of charge neutrality, this expression can be computed from Ewald potential Vew 
and the Madelung constant Vm, 

^ = E Q^Qj ^Ew(r^,) + ^ Q- Fm , (A2) 

i>j i 

where Vm = 5 limr-^o[^w('')~^c('')]- The Ewald potential is split in a real space and a Fourier part \§\ by subtracting 
the screening potential T4(r) of a Gaussian charge distribution. Under the assumption of charge neutrality Qi = 0) 
or in the presence of a neutralizing background, the Ewald potential becomes, 

VEw(r) = ^ [K - (r + L) + ^ K.(k)e»''-- + Cy , (A3) 

where the constant Cv, 

Cv^-^jd^v[V,-Vs]{v), (A4) 

represents contributions from a neutralizing background. For neutral systems, all Cy terms cancel. 

The total potential energy can be computed more efficiently by introducing the Fourier transform of the charge 
density, /5(k) = Qj e~'^^^i . Eg. lAll then reads. 



where the constant Vlmagc has been introduced. 



- VA (r + L) + \cv |p(0)p + \Y1 lp(k)l 



2 



^ ^ Qi Mmago , (-^^5) 



2Mmage = 2Vm " - ^ V;(k) = Y^[V,- V,] (L) - T4[0] . (A6) 

k^O Lt^O 

In contrast to Vm , Mmage depends on the choice of Vg ■ 

The same approach can be applied to construct the diagonal Ewald action where Vc (r) is then replaced by tVc (r) -|- 
Au(r, r; r). In many-body Monte Carlo simulations, an efhcient computation of the Ewald action, Eq. [211 is crucial 
This can be achieved by using the optimized Ewald break-up technique developed by Natoli et al. 



0, 



n |k|<fcc 
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where one uses only one real-space image and a variable number of Fourier vectors depending on the required accuracy. 
As basis set /„, we use locally piecewise quintic Hermite interpolants as suggested in However, this basis set is 
only used to perform the break-up. 

The coefficients a„ and are chosen to minimize the mean squared deviation between the UEw(r) and the fit 
function. While Natoli et al. proposed using a sum over many Fourier components for this step, we reached a higher 
accuracy by performing this calculation in real space, which leads to the usual set of n linear equations, 



= UEW.n - ^ {tEw(k)/„(k) - ^a„j 



k|<fcc 

with 



fm.n ^ ] /m(k)/n(k) 



k|<fcc 



(A8) 



UEW,n^ Jd^ruEw{^)fm{r) , /„,m = ^y"c^^r/n(r)/m(r) , (A9) 

and /„ and uew being the corresponding Fourier transforms of /„ and wew- The Fourier coefficients j/k are given by, 

2/k = ■iiEw(k) - ^ an /„(k) . (AlO) 

n 

The background term, Cu^ is treated as a A; = component. In the final step, the real space part, 

PF(|r|)=^a„/„(|r|) , (All) 

n 

is conveniently stored on a radial grid and interpolated during PIMC simulations. The same procedure is applied to 
the T derivative of the action. 
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TABLE I: Density matrix at = 8 for an isolated proton-electron pair. The columns correspond to diagonal potential 
action, the first off-diagonal expansion coefficient (see Eq. [14] with n = 1), the diagonal r derivative, and its first expansion 
coefficient. 



r 


u{r,0) 


Air) 


u'{r, 0) 




A' 


(r) 








-9 


052629e- 


-01 


-1.0656e-H00 


-3.701269e-H00 


4.6908e-H00 





1 


-7 


094574e- 


-01 


-7 


8668e- 


-01 


-3.664297e-H00 


2.2366e-h00 





2 


-5 


375241e- 


-01 


-5 


4357e- 


-01 


-3.435818e-H00 


4 


3544e- 


-01 





3 


-4 


052122e- 


-01 


-3 


5276e- 


-01 


-2.988023e-H00 


-5 


3328e- 


-01 





4 


-3 


130127e- 


-01 


-2 


1872e- 


-01 


-2.468529e-H00 


-8 


1126e- 


-01 





5 


-2 


512131e- 


-01 


-1 


3263e- 


-01 


-2.019987e-H00 


-7 


2479e- 


-01 





6 


-2 


090487e- 


-01 


-8 


0229e- 


-02 


-1.683282e-H00 


-5 


3529e- 


-01 





7 


-1 


789567e- 


-01 


-4 


6185e- 


-02 


-1.438155e-H00 


-3 


5572e- 


-01 





8 


-1 


564692e- 


-01 


-2 


7139e- 


-02 


-1.255464e-H00 


-2 


1970e- 


-01 





9 


-1 


390225e- 


-01 


-1 


7645e- 


-02 


-1.114415e-H00 


-1 


4285e- 


-01 


1 





-1 


250863e- 


-01 


-1 


2256e- 


-02 


-1.002118e-h00 


-9 


8852e- 


-02 


1 


1 


-1 


136946e- 


-01 


-8 


9124e- 


-03 


-9.105148e-01 


-7 


1705e- 


-02 


1 


2 


-1 


042075e- 


-01 


-6 


7064e- 


-03 


-8.343272e-01 


-5 


3873e- 


-02 


1 


3 


-9 


618330e- 


-02 


-5 


1838e- 


-03 


-7.699461e-01 


-4 


1601e- 


-02 


1 


4 


-8 


930751e- 


-02 


-4 


0954e- 


-03 


-7.148140e-01 


-3 


2843e- 


-02 


1 


5 


-8 


334981e- 


-02 


-3 


2949e- 


-03 


-6.670655e-01 


-2 


6411e- 


-02 


1 


6 


-7 


813769e- 


-02 


-2 


6922e- 


-03 


-6.253068e-01 


-2 


1571e- 


-02 


1 


7 


-7 


353934e- 


-02 


-2 


2291e- 


-03 


-5.884752e-01 


-1 


7856e- 


-02 


1 


8 


-6 


945233e- 


-02 


-1 


8672e- 


-03 


-5.557459e-01 


-1 


4954e- 


-02 


1 


9 


-6 


579581e- 


-02 


-1 


5801e- 


-03 


-5.264687e-01 


-1 


2653e- 


-02 


2 





-6 


250516e- 


-02 


-1 


3493e- 


-03 


-5.001243e-01 


-1 


0803e- 


-02 


2 


1 


-5 


952805e- 


-02 


-1 


1616e- 


-03 


-4.762926e-01 


-9 


2992e- 


-03 


2 


2 


-5 


682170e- 


-02 


-1 


0073e- 


-03 


-4.546301e-01 


-8 


0633e- 


-03 


2 


3 


-5 


435076e- 


-02 


-8 


7928e- 


-04 


-4.348534e-01 


-7 


0381e- 


-03 


2 


4 


-5 


208581e- 


-02 


-7 


7215e- 


-04 


-4.167263e-01 


-6 


1803e- 


-03 


2 


5 


-5 


000210e- 


-02 


-6 


8181e- 


-04 


-4.000506e-01 


-5 


4569e- 


-03 


2 


6 


-4 


807872e- 


-02 


-6 


0507e- 


-04 


-3.846586e-01 


-4 


8426e- 


-03 


2 


7 


-4 


629784e- 


-02 


-5 


3947e- 


-04 


-3.704075e-01 


-4 


3174e- 


-03 


2 


8 


-4 


464419e- 


-02 


-4 


8304e- 


-04 


-3.571749e-01 


-3 


8657e- 


-03 


2 


9 


-4 


310460e- 


-02 


-4 


3425e- 


-04 


-3.448554e-01 


-3 


4752e- 


-03 


3 





-4 


166768e- 


-02 


-3 


9181e- 


-04 


-3.333576e-01 


-3 


1355e- 


-03 



TABLE II: Pair density matrix for r ^ = 8 for isolated pair of electrons in the format of Tab. IH 



r 


u{r,0) 


A{r) 


u'(r,0) 


A'{r) 








6 


176418e- 


-01 


4 


4659e- 


-01 


2.435490e-(-00 


-1.6593e-|-00 





1 


5 


191880e- 


-01 


3 


3730e- 


-01 


2.424522e-|-00 


-9 


0614e-01 





2 


4 


290318e- 


-01 


2 


4933e- 


-01 


2.358555e-|-00 


-3 


6651e-01 





3 


3 


522642e- 


-01 


1 


8106e- 


-01 


2.218501e-f00 


-2 


2177e-02 





4 


2 


903868e- 


-01 


1 


2962e-01 


2.021059e-|-00 


1 


6609e-01 





5 


2 


422795e- 


-01 


9 


1878e- 


-02 


1.799564e-|-00 


2 


4375e-01 





6 


2 


0551646- 


-01 


6 


4773e- 


-02 


1.584616e-|-00 


2 


5324e-01 





7 


1 


774303e- 


-01 


4 


5605e- 


-02 


1.394356e-|-00 


2 


2765e-01 





8 


1 


557096e- 


-01 


3 


2174e-02 


1.234499e-|-00 


1 


8896^01 





9 


1 


385903e- 


-01 


2 


2809e- 


-02 


1.103098e-|-00 


1 


4953e-01 


1 





1 


248157e- 


-01 


1 


5221e- 


-02 


9.952834e-01 


1 


1214e-01 


1 


1 


1 


135143e- 


-01 


1 


0485e- 


-02 


9.060434e-01 


8 


1349e-02 


1 


2 


1 


04082 le- 


-01 


7 


6096e- 


-03 


8.312447e-01 


6 


0020e-02 


1 


3 


9 


eOgSWe- 


-02 


5 


7381e- 


-03 


7.677422e-01 


4 


5517e-02 


1 


4 


8 


924103e- 


-02 


4 


4539e- 


-03 


7.131934e-01 


3 


5419f:^02 


1 


5 


8 


329967e- 


-02 


3 


5366e- 


-03 


6.658461e-01 


2 


8165e-02 


1 


6 


7 


809915e- 


-02 


2 


8606e- 


-03 


6.243711e-01 


2 


2803e-02 


1 


7 


7 


350923e- 


-02 


2 


3497e- 


-03 


5.877450e-01 


1 


8744e-02 


1 


8 


6 


942845e- 


-02 


1 


9556e- 


-03 


5.551675e-01 


1 


5608e-02 


1 


9 


6 


577663e- 


-02 


1 


6462e- 


-03 


5.260046e-01 


1 


3144e-02 


2 





6 


248957e- 


-02 


1 


3996e- 


-03 


4.997475e-01 


1 


1178e-02 


2 


1 


5 


951525e- 


-02 


1 


2005e- 


-03 


4.759834e-01 


9 


5902e-03 


2 


2 


5 


681109e- 


-02 


1 


0377e- 


-03 


4.543740e-01 


8 


2918e-03 


2 


3 


5 


434190e- 


-02 


9 


0340e- 


-04 


4.346395e-01 


7 


2197e-03 


2 


4 


5 


207834e- 


-02 


7 


9148e- 


-04 


4.165462e-01 


6 


3261e-03 


2 


5 


4 


999577e- 


-02 


6 


9745e- 


-04 


3.998978e-01 


5 


5751e-03 


2 


6 


4 


807331e- 


-02 


6 


1784e- 


-04 


3.845282e-01 


4 


9393e-03 


2 


7 


4 


629319e- 


-02 


5 


4998e- 


-04 


3.702955e-01 


4 


3971e-03 


2 


8 


4 


464018e- 


-02 


4 


9176e- 


-04 


3.570782e-01 


3 


9319e-03 


2 


9 


4 


310111e- 


-02 


4 


4153e- 


-04 


3.447714e-01 


3 


5305e-03 


3 





4 


166463e- 


-02 


3 


9794e- 


-04 


3.332844e-01 


3 


1821e-03 
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TABLE III: Real space part W of the action for proton-electron and electron-electron pair density matrices computed for 
= 8 using an optimized Ewald break-up for a cubic simulation cell of size L = 5. The fc-space part is given in Tab. II VI 



r 




Wpe 


V } 






H4o 


(r) 


W (r) 


U 


U 


-3 


Q48Q72P- 


-01 


-2 4fi2finsp-+-on 


3 


779454P- 


-01 


1 Ti^Li*^ 1 Li Li 


n 
U 


i 


-2 


350750?- 


-01 


-2 258n2qp-i-on 


2 


9192fi0p- 


-01 


1 19fi514p-l-00 

i . -L i_iLIt_i -L^Cn^LiLi 


n 
U 


o 
z 


-1 


1 60451 p- 


-01 


-1 8525q3p-l-0n 


2 


1 90901 p- 

J- 0\J >J\J -Lc 


-01 


1 083259p-l-nn 

A • LiLJiJii^Jtycn^LiLi 


U 


Q 
O 


-3 


7SQ284P- 


-02 


-1 37f)n5fip-i-on 


1 


fi0849fip- 

LJ L/ (J t: Li C 


-01 


9 211 3n9p-ni 


n 
U 


/I 

4 


g 


20235QP- 


-03 


-9 fi33273p-01 


1 


1 fi3420p- 


-01 


7 390493p-ni 


n 
U 


cr 



3 


255803P- 

iJ O VJ <J c 


-02 


-fi fiOl 1 1 Ip-01 


g 


329489P- 


-02 


5 fi47894p-ni 


U 


/:? 
U 


4 


338204P- 


-02 


-4 51 1075P-01 


5 


925357P- 


-02 


4 1 557fi5p-ni 


U 


( 


4 


filfi727p- 


-02 


-3 0820fi7p-01 


4 


188fil4p- 

_L O Li -L t: ^3 


-02 


2 97291 3p-ni 

^ > C 1 ij<y i tJC Li -L 


U 


o 



4 


420Q1 7p- 


-02 


-9 09fifi09p-01 


2 


93974fip- 


-02 


9 nsi S89p-ni 


U 


9 


3 


Q4800SP- 


-02 


-1 41 31 32p-01 


2 


041 829p- 


-02 


1 432n97p-fll 


i 


U 


3 


3387f!fip- 


-02 


-9 395882P-02 


1 


400554P- 

t: Li Li t_i t: C 


-02 


9 fi83finfip-n2 

.LiOtJLiLiLiC LiiJ 


i 


i 


2 


70345SP- 


-02 


-fi 140030P-02 


9 


4fifil 5fip- 


-03 


fi 429422p-fl9 


i 


Z 


2 


11037Qe- 


-02 


-3.929014e-02 


6 


284010e- 


-03 


4 181878e-02 


1 


Q 

o 


1 


5875f!Qp- 


-02 


-2 4535fi8p-02 


4 


085855P- 

Li C CJ !_i t_i c 


-03 


2 fi57455p-n2 


1 
i 


/I 




1 497SQp- 


-02 


-1 489252P-02 


2 


581 234p- 


-03 


1 fi441 21 p-fl9 


1 
i 


r 



7 


866995e- 


-03 


-8 749983e-03 


1 


587436e- 


-03 


9 870341e-03 


I 


5 


5 


231 1 15p- 


-03 


-4 948775P-03 


9 


4270fi2p- 

1 LiLi^Ci 


-04 


5 723989p-n3 

t_i . 1 ^tJtJUiJy^ LitJ 


1 


7 


3 


356079e- 


-03 


-2.677517e-03 


5 


456568e- 


-04 


3.186790e-03 


1 


8 


2 


013496e- 


-03 


-1.378287e-03 


3 


001534e- 


-04 


1.692584e-03 


1 


9 


1 


098643e- 


-03 


-6.634200e-04 


1 


574300e- 


-04 


8.477546e-04 


2 





5 


739410e- 


-04 


-2.995616e-04 


7 


289096e- 


-05 


3.991008e-04 


2 


1 


3 


066066e- 


-04 


-1.223001e-04 


2 


820100e- 


-05 


1.695370e-04 


2 


2 


1 


462296e- 


-04 


-4.510938e-05 


1 


383209e- 


-05 


6.587290e-05 


2 


3 


3 


113736e- 


-05 


-1.503693e-05 


3 


426564e- 


-06 


1.986137e-05 


2 


4 


-1 


417605e- 


-05 


-4.564337e-06 


6 


935208e- 


-07 


6.273595e-06 


2 


5 
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TABLE IV: K-space action for proton-electron and electron-electron pair density matrices computed for = 8 using an 
optimized Ewald break-up with 20 k-shells for a cubic simulation cell of size L — 5. The corresponding real space parts are 
listed in Tab. IIIll The table also includes the Madelung and the background constant from Egs. 1251 and IA7I 



n^, k = 2%n/L 






"^fc ,CC 






3.5466e- 


-02 


2.8372e- 


-01 


-3.5467e-02 


2.8376e-01 


Cu 


-3.0502e- 


-03 


2.0762e- 


-02 


-2.9417e-03 


-2.2041e-02 


1 


-9.8752e- 


-03 


-4.5714e- 


-02 


5.4742e-03 


4.4979e-02 


2 


-5.0783e- 


-03 


-1.6276e- 


-02 


1.8624e-03 


1.5781e-02 


3 


-3.1684e- 


-03 


-7.6607e- 


-03 


8.3439e-04 


7.3217e-03 


4 


-2.0962e- 


-03 


-4.0198e- 


-03 


4.1473e-04 


3.7862e-03 


5 


-1.4182e- 


-03 


-2.2280e- 


-03 


2.1645e-04 


2.0676e-03 


6 


-9.6702e- 


-04 


-1.2727e- 


-03 


1.1557e-04 


1.1634f:^03 


8 


-4.4850e- 


-04 


-4.3268e- 


-04 


3.3300e-05 


3.8330e-04 


9 


-3.0283e- 


-04 


-2.5364e- 


-04 


1.7587e-05 


2.2113e-04 


10 


-2.0262e- 


-04 


-1.4824e- 


-04 


9.0566e-06 


1.2710e-04 


11 


-1.3405e- 


-04 


-8.5951e- 


-05 


4.4823e-06 


7.2484e-05 


12 


-8.7505e- 


-05 


-4.9194e- 


-05 


2.0841e-Q6 


4.0853e-05 


13 


-5.6226e- 


-05 


-2.7704e- 


-05 


8.7102e-07 


2.2662e-05 


14 


-3.5469e- 


-05 


-1.5293e- 


-05 


2.9069e-07 


1.2311(:^05 


16 


-1.3175e- 


-05 


-4.2131e- 


-06 


-5.3452e-08 


3.3414e-06 


17 


-7.6876e- 


-06 


-2.0560e- 


-06 


-7.1360e-08 


1.6470e-06 


18 


-4.3189e- 


-06 


-9.4808e- 


-07 


-6.0744e-08 


7.6995e-07 


19 


-2.3129e- 


-06 


-4.0035e- 


-07 


-4.2831e-08 


3.3463e-07 


20 


-1.1631e- 


-06 


-1.3415e- 


-07 


-2.6457e-08 


1.3167e-07 


21 


-5.3551e- 


-07 


-1.5653e- 


-08 


-1.4397e-08 


4.4726f:^08 


22 


-2.1498e- 


-07 


1.9113e- 


-08 


-6.7088e-09 


1.1253e-08 



